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ABSTRACT 

Context. Hot molecular cores (HMCs) are intermediate stages of high-mass star formation and are also known for their rich chemical 
reservoirs and emission line spectra at (sub-)mm wavebands. Complex organic molecules (COMs) such as methanol (CH 3 OH), ethanol 
(C 2 H 5 OH), dimethyl ether (CH 3 OCH 3 ) and methyl formate (HCOOCH 3 ) produce most of these observed lines. The observed spectral 
feature of HMCs such as total number of emission lines and associated line intensities are also found to vary with evolutionary stages. 
Aims. We aim to investigate the spectral evolution of these COMs to explore the initial evolutionary stages of high-mass star formation 
including HMCs. 

Methods. We developed various 3D models for HMCs guided by the evolutionary scenarios proposed by recent empirical and mod¬ 
eling studies. We then investigated the spatio-temporal variation of temperature and molecular abundances in HMCs by consistently 
coupling gas-grain chemical evolution with radiative transfer calculations. We explored the effects of varying physical conditions 
on molecular abundances including density distribution and luminosity evolution of the central protostar(s) among other parameters. 
Finally, we simulated the synthetic spectra for these models at different evolutionary timescales to compare with observations. 
Results. Temperature has a profound effect on the formation of COMs through the depletion and diffusion on grain surface to 
desorption and further gas-phase processing. The time-dependent temperature structure of the hot core models provides a realistic 
framework for investigating the spatial variation of ice mantle evaporation as a function of evolutionary timescales. We find that 
a slightly higher value (15K) than the canonical dark cloud temperature (lOK) provides a more productive environment for COM 
formation on grain surface. With increasing protostellar luminosity, the water ice evaporation font ('-lOOK) expands and the spatial 
distribution of gas phase abundances of these COMs also spreads out. We calculated the temporal variation of the radial profiles of 
these COMs for different hot core models. These profiles resemble the so-called jump profiles with relative abundances higher than 
10““* within the evaporation font will furthermore be useful to model the observed spectra of hot cores. We present the simulated 
spectra of these COMs for different hot core models at various evolutionary timescales. A qualitative comparison of the simulated and 
observed spectra suggests that these self-consistent hot core models can reproduce the notable trends in hot core spectral variation 
within the typical hot core timescales of 10^ year. These models predict that the spatial distribution of various emission line maps will 
also expand with evolutionary time; this feature can be used to constrain the relative desorption energies of the molecules that mainly 
form on the grain surface and return to the gas phase via thermal desorption. The detailed modeling of the thermal shucture of hot 
cores with similar masses along with the characterization of the desorption energies of different molecules can be used to constrain 
the luminosity evolution of the central protostars. The model predictions can be compared with high resolution observation that can 
probe scales of a few thousand AU in high-mass star forming regions such as from Atacama Farge Millimeter/submillimeter Array 
(ALMA). We used a spectral fitting method to analyze the simulated spectra and find that it significantly underestimates some of the 
physical parameters such as temperature. The coupling of chemical evolution with radiative transfer models will be particularly useful 
to decipher the physical structure of hot cores and also to constrain the initial evolutionary stages of high-mass star formation. 
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1. Introduction 

Hot molecular cores (HMCs), also known as hot cores, are asso¬ 
ciated with the early evolutionary stages of high-mass star for¬ 
mation. Our current understanding of high-mass star formation 
is incomplete, but a plausible evolutionary scenario, emerging 
from the recent observational studies, is as follows; Infrared dark 
clouds (IRDCs) ^ hot molecular cores —» hyper/ultra-compact 
Hit regions —> Hii region surrounding the ionizing high-mass 
stars (see Beuther et al.pOOTl ). The intermediate phases between 
IRDCs and hypercompact Hit regions are known as hot cores; 


at this phase the central protostars heat up their surrounding 
gaseous and dusty materials to high temperatures (> 100 K) 


rchoudhury @ mpe. mpg. de 


without ionizing it (due to lack of UV photons). Kurtz et al. 
( |2000) l characterized HMCs as a relatively compact phase of 
high-mass star formation (< 0.1 pc) with a typical density and 
temperature in the range of > 10^ cm~^ and > 100 K. The evapo¬ 
ration of water and organic rich ice mantles leads to the detection 
of a wide variety of complex organic molecules (COMs) such as 
CH3O H, CH3OCH3, HCOOCH3 and C2H5OH in a number of 
HMCs ( |Herbst & van Dishoeck|2009l l. 

The chemical richness and diversity of HMCs has received 
significant attention in recent times and has also triggered a num- 
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ber of studies that modeled the relevant chemical processes re¬ 
sponsible for building up the molecular reservoir in HMCs (see 


Herbst & van Dishoeck 2009 for a detailed review of existing 
models). Among these, gas-grain models (which include both 
gas phase and grain surface chemical reactions) emerged as the 
most promising ones to explain the abundance of COMs in hot 
cores ( |Charnley et al.|[T992| l. These studies suggested that the 
chemical richness of HMCs originates from the evaporation and 
subsequent gas phase processing of the molecules that became 
locked on icy mantles during the initial collapse phase of high- 
mass protostar formation. Recent hot core models (Garrod et al. 
2008 |Garrod|2013a| l successfully reproduced the abundances of 
a number of COMs inferred from the observational studies of hot 
cores. One of the crucial parameters of these models is the grad¬ 
ual warm-up of hot cores to a temperature of 100 K and higher 
to facilitate the formation of COMs on grain surface and even¬ 
tually their evaporation to the gas phase. Recent observations 
indicate that COMs are also present in cold dark clouds with a 
typical temperature of 10 K ( Cernicharo et al.||2012[ |Bacmann| 
et al.|20T2| . Vasyunin & Herbst (20131 suggested some alter¬ 
native mechanisms that do not require high temperature (e.g., 
reactive desorption) for COM formation in these environments. 
However, these processes are probably less important in the hot 
core phase. 

The salient observational feature of HMCs are numerous 
emission lines at (sub-)mm wavebands, mostly originating from 
the rotational transitions of COMs. Dedicated line surveys (e.g.. 


Schilke et al.||1007| 

20011 [Zemickel et al.||2012[ [Crockett et al. 

20141 Neill et al.||20141 see also Table. 2 of IHerbst & van 

Dishoeckj 2009 

1 are very useful to explore the inventory of 


COMs in hot cores. Most of these studies used local thermo¬ 
dynamic equilibrium (LTE) fitting of spectra, for example, rota¬ 
tional diagram analysis ( [Goldsmith & Langer||1999) , XCLASS 
(Zemickel et al.|2012|), or similar techniques for estimating the 


physical parameters and chemical abundances. However, these 
techniques, while useful to estimate the average physio-chemical 
scenario of the hot cores, are not suitable to explore the de¬ 
tailed structures (Herbst & van Dishoeck [2009 |l. Investigations 


of the chemical evolution of selected molecules using a num¬ 
ber of hot core sources are also carried out in selected studies 
( Bisschop et al.||2007 Beuther et al.||2009 Gemer et al.|2014 1 . 
It has been found that the intensity as well as the spatial distri¬ 
bution of emission lines of various molecules also varies with 
evolutionary stages and from one hot core to another. These fea¬ 
tures can be attributed to the varying physio-chemical condition 
during the initial evolutionary stages of high-mass star forma¬ 
tion. Although the results are promising, high angular resolution 
spectral line observations of the initial stages of high-mass star 
formation are still rare ( jBrouillet et al.|20T3| ). 

With the availability of the Atacama Large Millime¬ 
ter/submillimeter Array (ALMAQ it is now possible to ex¬ 
plore the chemical segregation of high-mass star forming regions 
down to a spatial resolution of a few hundred AU in relatively 
short time. Observations of many sources at various evolution¬ 
ary stages will also provide the opportunity to explore the chem¬ 
ical diversity. In short, these observed chemical probes can be 
used as powerful tools to decipher the internal structure and to 
constrain the evolutionary stages of hot cores. So far, only few 
astrochemical studies (e.g., Doty et al.|2002 Qberg et al.|2013 


Gerner et al.|2014 1 incorporated a radial variation of density and 
temperature. However, a self-consistent approach to explore the 
temporal evolution of physical parameters associated with high- 


* http://www.almaobservatory.org/ 
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Fig. 1: Flowchart of the modeling framework. 


mass star formation is not considered in these models. Moreover, 
there are no dedicated studies that explore the trends of tempo¬ 
ral variation in simulated spectra which arise due to the vary¬ 
ing physical and chemical conditions of evolving hot cores. To 
bridge this gap, we investigated the effects of the spatio-temporal 
variation of physical conditions on the chemical evolution of hot 
cores in a self-consistent way. In this work, we investigate the 
chemical evolution of selected oxygen-bearing COMs such as, 
CH 3 OH, CjHjOH, HCOOCH 3 and CH 3 OCH 3 and simulate the 
spectral data cubes at (sub-)mm wavebands for various hot core 
models at different evolutionary stages. 

We describe the model set-up, initial physical conditions 
and chemical modeling in Sect. Spatio-temporal evolution 
of molecular abundances of selected COMs are summarized in 
Sect.[^ The simulated spectra and a qualitative comparison with 
the observed spectra are presented in Sect.|^ In Sect.|^ we sum¬ 
marize this study and discuss the possibilities of using these sim¬ 
ulated spectra along with the observations to explore the evolu¬ 
tion of hot cores, and finally we present our conclusion in Sect.|^ 


2. Method 

Most of the high-mass star forming regions are distant and heav¬ 
ily obscured. The physical structure, evolutionary sequence and 
timescales of high-mass star formation are thus only poorly un¬ 
derstood. Molecular emission lines provide good diagnostics to 
understand the internal structure of these environments ( jRolffsj 
jet al.|20lT] l. In the case of low-mass star formation, several stud¬ 
ies incorporated the results of hydrodynamic simulations into 
chemical models to investigate the role of evolving physical con¬ 
ditions on the chemical composition ( jAikawa et al. [200^ [Fumy a] 
jet al.|201^ . Similarly, the ideal way to model hot core chemistry 
would be to couple the chemical modeling with hydrodynami- 
cal simulations of high-mass star formation. However, models 
of high-mass protostellar evolution along with the surrounding 
gaseous and dusty materials are rare, with only a handful of 
available solutions ( jKuiper & Yorke||2013| l. In order to inves¬ 
tigate the coupling of physical conditions and molecular abun¬ 
dances, we set-up a simple physical model of hot cores guided by 
the current understanding of high-mass star formation and used 
this set-up to examine the spatio-temporal variation of chemical 
evolution of hot cores and finally generated synthetic spectra fol¬ 
lowing the flowchart shown in Fig. The flowchart is executed 
using a Python script (Schmiedeke et al., 2014, in preparation). 
We describe the physical parameters and chemical models in de¬ 
tail below. 
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2. 1. Chemical model 


We developed a rate-equation-based ID astrochemical code 
iSaptarsy) in Fortran90; it calculates the spatio-temporal evo¬ 
lution of various molecular species using gas phase reactions, 
gas-grain interactions, and surface chemistry (Choudhury et al., 
2014, in preparation). The code has its origin in the Astrochem 
code of Bergin et al. ( |1995[ ), and after several major modifica¬ 
tions, it has become an improved and new astrochemical code 
that can be easily integrated with other external codes (e.g., 
RADMC-3I^. The ordinary differential equation solver from 
Netlib library dvodpl^and MA2^ the advance solver of sparse 
systems of linear equations from the HSL library, are used to 
compute the chemical evolution. Gas-grain interaction via ac¬ 
cretion and desorption and surface reactions are implemented 
following the recipe of |Hasegawa et al.] ( |1992[ ) and |Hasegawa &| 
Herbst (|1993|l. The temporal evolution of grains such as grain 


coagulation and formation of the grain mantles etc. can affect 
the grain surface chemistry ( [Taquet et al.|2012] |Garrod|2013b| l, 
but these processes are beyond the scope of this study and are, 
therefore, not considered in this work. 

Saptarsy is suitable to explore the effects of spatio-temporal 
variation of physical parameters (e.g., temperature) on chemi¬ 
cal evolution. It has been benchmarked against the models de¬ 
scribed in [Semenov et al.| ( |2010| l. The gas-grain chemical net- 
worlj^from the Ohio State University Astrophysical Chemistry 
Group (~7840 reactions and ~840 species) is used in this study. 
The detailed description of the network and reaction rate calcu¬ 
lations of the associated chemical reactions can be found in IGar-l 
[rod et al.| ( |2008| l and |Semenov et al.] ( |2010[ ). We adopted standard 
grain parameters such as radius of 0.1 pm , density of 3 g cm“^ 
and a dust-to-gas mass ratio of 0.01. The initial abundance and 
desorption energy of the molecular species were adopted from 
[Garrod et al.|P0^08l l and |Garrod| ( |20 13 ^ . The ratio of diffusion- 
to-desorption energy is assumed to be 0.5. We used Saptarsy to 
obtain the spatio-temporal abundance evolution by solving the 
rate equations based on the overall molecular formation and de¬ 
struction terms. Chemical evolutionary models are simulated up 
to a time-scale of 1 x 10 ^ year using discrete time-steps sepa¬ 
rated by few hundreds to thousands year. Conservative values of 
relative and absolute tolerances of 10 “^ and 10 “^°, respectively, 
were used throughout the calculation. To avoid complexity, di¬ 
rect photo-dissociation and -ionization reactions are not con¬ 
sidered in this work since hot core phases were not dominated 
by UV radiation. However, the effects of UV photons become 
prominent at the final phases of hot core evolution that is during 
the transformation to ultracompact Hit regions which require a 
detailed modeling of the radiation field inside hot cores; this is 
the focus of a separate study (Stephan et al., 2014, in prepara¬ 
tion). 


2.2. Density distribution 

In a complex environment such as a hot core, the physical condi¬ 
tions (e.g., density, temperature) not only vary with time, but also 
with the distance from the central protostars. Recent (sub-)mm 
interferometric studies also indicated that internal rotational mo¬ 
tion, outflows, etc. can be associated with HMCs ( [Beltran et al.| 
[201 1[ [Furuya et al. 12011] ). Several studies investigated the den¬ 
sity distribution of these objects using dust continuum emission 

^ http://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/ 

^ http://www.netlib.org/ode/ 
http://www.hsl.rl.ac.uk/ 

^ http://www.physics.ohio-state.edu/~eric/research.html 


Table 1: List of abbreviations used in the nomenclature of vari¬ 
ous models 


Description/parameter 

Value 

Abbrv. 

Cold core models 

_ 

cc 

Hot molecular core models 

- 

hmc 

Plummer radius 

3000 AU 

r3 


3500 AU 

r3.5 


4000 AU 

r4 


4500 AU 

r4.5 

Constant temperature 

10 K 

tlO 

at cold core phase 

15 K 

tl5 

Luminosity evolution 

L^ 

12 

during warm-up phase 

L4 

14 


v-^ 

11 

Cosmic-ray ionization 

1.37 X 10-1'^ s-i 

crl7 

rate 

1.37 X 10-'® s-' 

crl 6 


(van der Tak et al. 2000 Beuther et al.|2002) l and found a power- 
law density distribution on larger scales; recent studies using in¬ 
terferometric observations suggested a Plummer-type distribu¬ 
tion at smaller scales ( Rolffs et al.|2011|l . Plummer-type profiles 
can be represented as n(r)= nc(l-Hc(r/rp)^)'>', where c is a constant 
and n^ and r^ represent the central density and Plummer radius 
(this parameter sets the size of the core at which n,, drops to half 
of its value) and the exponent (y) sets the density distribution at 
radii larger than the Plummer radius. We set up a spherical hot 
core model that contains a high-mass protostar at the center em¬ 
bedded in a dense core of gas and dust. We adopted a Plummer- 
type distribution for the hot core models with a peak density of 
2 X 10^ cm“^, c=0.31951 and y=-5/2 throughout this work. The 
density distribution profile does not change with time, but we 
have considered various types of density profiles by varying the 
Plummer radius (see Table. T]|. Assuming the standard dust-to¬ 
gas mass ratio, the mass of the model core with Plummer radius 
~3000 AU is estimated to be around 45 Mq, indicating that it 
barely has the potential to form a high-mass star (> IOMq). 


2.3. Protostellar evolutionary models 

The central protostar is the only source of heating in our model 
set-up and thus the temperature of hot cores depend on the lu¬ 
minosity of the protostar. We used the protostellar evolution¬ 
ary models of [Hosokawa & Omukai[ ( [2009| l (spheric al accretion 
model with a constant accretion rate ofTO^ Mgyr ' ([Rolffs et al. 


201 1 ||) to investigate the spatio-temporal evolution of the tern 


perature in our hot core models. The temperature evolution from 
IRDC phases (10-15 K) to the formation of protostars is not 


explicitly covered in this model. Garrod et al. (20081 explored 


the effect of warm-up timescales (the required time to reach the 
characteristic temperature of 100 K) and found that with a longer 
warm-up timescale the chemical complexity also increases. To 
imitate the luminosity evolution during the warm-up phase, we 
interpolated the luminosity over a fixed time interval (5x10^ to 
7.5 X 10"^ year) using different power-law indices of 2, 4, and 7.5 
(see Fig.j^a)). The starting point of the luminosity interpolation 
(5 X 10"^ year) was chosen in such a way that the temperature at 
the immediate vicinity of the protostars is around 15 K; the lumi¬ 
nosity at the end point is adopted from the models of |Hosokawa[ 
& Omukai so that the temperature reaches up to 100 K at the 


center. We then used these luminosity profiles to calculate the 
temperature. The power-law index 7.5 was chosen to reproduce 
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a similar temperature evolution as used in |Garrod et al.| ( [2008| l. 
We restricted the highest value of the protostellar luminosity to 
1 X 10^ Lq following the result of recent observational studies of 
hot cores ( jBeuther et al.|2009] l. 

2.4. Nomenclature of the models 

We developed various models by varying the physical param¬ 
eters to explore the spatio-temporal variation of the chemical 
evolution. Instead of commonly used tags (e.g., modeha), we 
use different combinations of abbreviations to refer to a partic¬ 
ular model. We summarize the physical parameters and associ¬ 
ated abbreviations in Table For example, a hot core model 
tagged as crl6-r3-l2 refers to a model where the Plummer ra¬ 
dius is 3000 AU (r3), the power-law index of the protostel¬ 
lar luminosity evolution is 2 (12), and the cosmic-ray ioniza¬ 
tion rate is 1.37 x 10“'® s“*. Similarly, the cr-17 models refers 
to all the hot core models with a cosmic-ray ionization rate of 
1.37 X 10“'^ s-'. 


2.5. Simulation of spectral data cubes 


The radiative transfer code RADMC-3D was used to simulate 
the spectral data cubes of different COMs at various wavebands 
and evolutionary stages; the radius adopted for the hot core mod¬ 
els was 0.1 pc (~20600 AU). We started with a density distribu¬ 
tion and a luminosity evolution prohle of the central protostar 
and used RADMC-3D to calculate the spatio-temporal evolu¬ 
tion of the temperamre. In the next step, we extracted the ID 
density and temperature distribution profiles from the 3D hot 
core models. These density-weighted profiles along with the 
initial abundance of selected molecules served as the input for 
Saptarsy, which calculates the spatio-temporal chemical evo¬ 
lution in hot cores. The ID abundance distribution prohles of 
respective molecules were then used as input abundance pro¬ 
hles in RADMC-3D. Finally, simulated spectra spanning over a 
frequency range at desired evolutionary timescales were gener¬ 
ated assuming local thermodynamic equilibrium (LTE). As sug¬ 
gested by |Herbst & van Dishoeck ( |2009 1, the ideal way to gen¬ 
erate simulated spectra would be to use full non-LTE statisti¬ 
cal equilibrium excitation of the molecules. However, collisional 
rate coefficients required for non-LTE modelings are available 
only for a few molecules. Moreover, the high densities and ra¬ 
diative coupling with far-infrared radiation helds make LTE a 
good approximation. Therefore, for this study we only consid¬ 
ered LTE excitation to generate comparable simulated spectra 
for different COMs. To facilitate the comparison with observa¬ 
tions, simulated spectra were convolved with different telescope 
beam sizes using MIRIAI^ also taking into account the distance 
to the source. The hot core models were placed at a hducial dis¬ 
tance of 2 kpc. Molecular data hies of various COMs were re¬ 
trieved fro m the Cologne Database fo r Molecular Spect roscopy 
(CDMSfllMuller et al.||2005| |200l| l, t he JPL Catalog (|Pickett 
|et al.|19M| l, and the VAMDC database ( |Rixon et al.|201 l^ T 


3. Results 

3.1. Molecular abundance on grains 

The abundance of various species on grains that is the ice abun¬ 
dances, are important parameters that inhuence the chemical 

® http://bima.astro.umd.edu/miriad/ 

’ http://cdms.phl.uni-koeln.de/cdms/portal/ 
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evolution of hot cores. At low temperature (10-15 K), molecules 
accrete onto grains and thus become locked on the grain surface; 
this is also known as freeze-out or depletion of molecules. We 
represent the grain surface species with s- meaning that s-HjO 
represents water on grain surface. A reduced abundance of some 
molecules at the center of the cores, that is in the high density 
regions, have been observed in a number of starless cores indi¬ 
cating the clear presence of molecular ice formation at this stage 
( jBergin & Tafalla||2007| l. The accretion rate of molecules onto 
grain surfaces depends on the sticking coefficient of molecules 
to the grain, the grain surface area, the velocity of the molecules 
in gas phase, and the number density of dust grains. Another 
critical factor that regulates the ice formation is the dust temper¬ 
ature at the early pre-stellar stages. Most of the hot core models 
so far considered the initial temperature as 10 K. However, re¬ 
cent large scale surveys of high-mass star forming regions indi¬ 


cate an average temperature of 15 K for starless cores (Wienen 
et al.|201^ Hoq et al.|2013|l. The slight increase of temperature 


from 10 to 15 K increases the mobility of molecules on grains 
(reactive radicals such as OH, CHj) and thus provides a more 
productive environment to form complex molecules. 

We therefore explored models with temperatures of 10 K and 
15 K during the prestellar phase. Another factor that influences 
the abundance of molecules frozen to grain surfaces is the inter¬ 
action with cosmic rays that can dissociate and even evaporate 
ices. There is a debate about the cosmic-ray rate in the lite rature 
based on the observations of H 3 (Indriolo & McCall 2012 1 . Mo¬ 
tivated by these studies, we also explored the ettect of cosmic- 
ray ionization rate wi th two distinct values of. 1.37 x 10“'® s“' 
and 1.37 x 10“'^ s ' ( Padovani et al.|2009| . We ran four differ¬ 
ent models to calculate the ice abundances, cc-crl7-tl0, cc-crl6- 
tlO, cc-crl7-tl5, and cc-crl6-tl5. The timescale of the starless 
phase is uncertain, but we have some guidance from large scale 
survey s of high-mass star forming regions ( jTackenberg et ah] 
2012 1 ; thus we adopted a timescale of 5 x 10"' year for the pre- 


stellar phase. In Eig.j^we present the grain surface abundance of 
selected molecules at the pre-stellar stage with respect to S-H 2 O 
in percent. All these models show a similar s-HjO relative abun¬ 
dance 2 X 10“"' with respect to the total hydrogen number density 
(n(H-i- 2 H 2 )) at the center. We also tabulated the ice abundance 
of these molecules from the literamre in Table Comparing the 
observed and our simulated values, it appears that cold core mod¬ 
els cc-crl7-tl5 and cc-crl6-tl5 produce S-CH 3 OH abundances 
close to the observed values among the other molecules. Since 
S-CH 3 OH acts as the parent molecule of other COMs, we have 
used these two models for the subsequent calculations. 


3.2. Spatio-temporal variation of temperature 

Temperature plays an important role in chemical evolution as 
the rate of most of the chemical reactions strongly depend on it. 
Previous studies (e.g., |Viti & Williams|[l999[ |Garrod & Herbst| 
2006) 1 investigated the effect of temperature evolution assuming 
a power-law distribution in a constant density set-up. However 
no self-consistent treatment is used in these studies to calcu¬ 
late the spatial structure of temperature using any protostellar 
luminosity function. A significant improvement for temperature 
evolution is achieved by using a proper radiative transfer calcu¬ 
lation that is tied to an evolving protostellar luminosity model. 
As required by the chemical evolutionary model Saptarsy, tem¬ 
perature evolution of the hot core models for different lumi¬ 
nosity evolution are calculated at some discrete time-steps us¬ 
ing RADMC-3D (Pig|^. We estimated the heating and cooling 
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Fig. 2; Evolution of molecular abundance on grains with respect to S-H 2 O in percent at 5x10"^ year for different cold core models: 
(a) cc-crl7-tl0, (b) cc-crl6-tl0, (c) cc-crl7-tl5, and (d) cc-crl6-tl5. 


Table 2: 

Molecular abundances on grains with respect to 

5 -H 2 O in percent: observed and simulated values 

Molecules 

W33A^ 

NGC 7538 
IRS 92 

Sgr A*^ 

Garrod et al. 


Oberg et al. 


Garrod 

Model 

Model 

cc-crl6-tl5 


(2DD8) 

(2011) 

(2013) 

cc-crl7-tli 

CO 

8 

16 

<12 

19 

13 

57 

8.75 

7.9 

CO 2 

13 

22 

14 

4.1(-3) 

13 

18 

0.01 

0.005 

CH 4 

1.5 

2 

2 

22 

2 

1.9 

35.2 

36.6 

NH 3 

15 

13 

20-30 

25 

5 

18 

22.7 

23.1 

H 2 CO 

6 

4 

<3 

4.3 

<2 

1.6 

0.22 

9(-10) 

HCOOH 

7 

3 

3 

3.2(-6) 

- 

- 

0.21 

1.7(-10) 

CH 3 OH 

18 

5 

<4 

15 

4 

6.9 

1.03 

4.0 


References. (1) Gibb et al. 1 2000b 1 ; (2) see Gibb et al. 1 2000b 1 for original references. a(b) refers to ax 10* throughout the manuscript. 


timescales of dust grains explicitly to check whether dust tem¬ 
perature calculation at discrete time-steps would yield a realis¬ 
tic estimate of temperature evolution. To estimate the cooling 
timescales, we compared the absorption and emission processes 
of dust grains. A typical dust grain of radius 0.1 jim at 20 K 
would emit 1.14 x 10“'^ erg s“* assuming an efficiency fac¬ 
tor of 10 “^ (considering the opacities at (sub-)mm wavelengths). 
Similarly, considering our assumed density of grains (3 g cm“^), 
the thermal energy of the dust grain would be 2.5 x 10“^ erg. 
Therefore, the time-scale for thermalization of dust grains would 
be about the fraction of a year (private communication with C. 
P. Dullemond). As a result, any changes in temperature of the 


hot core model become nearly instantaneous compared to the 
changes in luminosity and thus ensure a nearly steady-state con¬ 
dition in the hot core models. Dust mass opacity coefficients 
were retrieved from |Ossenkopf & Henning] ( |1994| l. We also as¬ 
sumed that gas temperature is same as the dust temperature, 
which is a reasonable assumption considering the high density 
of hot cores. |Doty et ^ ( |2002[ ) made a detailed calculation of 
the gas temperature for their hot core model and found that the 
values are similar to the ones obtained using the assumption of 
Tgas-^dust- In Figs. J^b) and[^c), the spatio-temporal evolution 
of temperature for hot core models hmc-r3 are shown. At the 
center of the hot cores, the temperature remains constant (15 K) 
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Fig. 3; (a) Luminosity evolution profiles for different power-law indices, (b) Temperature evolution in the hot core model hmc-r3 
at different radii; ~200 AU (blue), ~2000 AU (green) and ~ 10000 AU (red). Temperature distribution obtained using different 
protostellar luminosity evolutions are shown with continuous (12), dashed (14), and dotted (17) lines, (c) Spatio-temporal evolution 
of the temperature in the hot core model hmc-r3. 


up to 5 X 10“* year (cold core phase), then it rises to the typi¬ 
cal hot core temperature (100 K) during the next 2.5 x 10"^ year 
(warm-up phase), and after that, it reaches up to 500 K at the 
end of the evolutionary time-scale that is 1 x 10 ^ year (hot core 
phase). It is evident from the figures that the duration of these 
phases varies as a function of distance from the central star. 


3.3. Spatio-temporal abundance variation 


The most abundant molecule that dominates the spectra at 
(sub-)mm wavebands is CH 3 OH. Following CH 3 OH, the other 
oxygen-bearing COMs that contribute most of the lines in (sub- 
)mm wavebands (480-1907 GHz) are CH 3 OCH 3 , HCOOCH 3 


and CjHjOH (jZernickel et al.|2012|l. These COMs are also re¬ 
ferred as first-generation molecules (Bisschop et al.|2007[[G^ 


[rod et al.|2008[ [Herbst & van Dishoeck|2009[ ); these molecules 

are mainly produced by reactions between primary ice species 
(also known as the zeroth-generation molecules). To inves¬ 
tigate the favorable physical conditions for the formation of 
first-generation COMs during the warm-up phase, we ran fur¬ 
ther mode ls using the initial physical conditions discussed in 


Sect. 3.1 We ran 24 hot core models in total up to 1 x 10^ 


year with different physical parameters such as with four dif¬ 
ferent static density distributions (r3, r3.5, r4, and r4.5) along 
with three different luminosity evolutions (12, 14, and 17) and 
two different cosmic-ray ionization rates (cr-16 and cr-17) (see 
Table [^. The dust temperature is one of the most dominant fac¬ 
tors that control the gas and grain surface abundances. Thanks 
to the spatio-temporal evolution of the temperature structure, 
our model set-up is particularly useful to investigate the tran¬ 
sition from grain to gas phase abundance of various molecules 
as a function of evolutionary stages and distances from the cen¬ 
tral star. The spatio-temporal variations of grain and gas-phase 
abundance of selected COMs for various hot core models are 
shown in Figs. A.l A.4 At the initial stages, the ice abun¬ 


dance of COMs dominate over the gas phase abundances but 
with increased temperature ices evaporate from the grain sur¬ 
faces and thus boost the gas phase abundances. With temporal 
evolution, the radius of the evaporation font (where the temper¬ 
ature is equal to the water ice evaporation temperature of 100 K) 
also increases (see Fig. ED. and consequently, the spatial dis¬ 
tribution of gas phase abundances also spreads out. However, 
at the outer region (>4.5 x 10^ AU), where the temperature re¬ 
mains lower than 100 K, ice abundance of COMs dominate the 
gas phase abundances. One of the important features of this work 
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is the derived temporal variation of the spatial distribution of 
COMs for solid and gaseous phases. These profiles will be par¬ 
ticularly useful for modeling the observed high resolution emis¬ 
sion line maps of these molecules. Gas phase abundance of se¬ 
lected COMs from the existing observational and modeling stud¬ 
ies are collected in Table|^ a similar compilation of the observed 
value of these molecules for a number of hot cores can be found 
in [Bisschop et al.| ( |2007| ). The gas phase abundance of different 
COMs at the center of various hot core models are summarized 
in Tables |^and|^ These abundances are comparable with the ob¬ 
served abundance of various COMs at least for the crl6 models. 
However, while comparing the observed and simulated values, 
one should also keep in mind that the quoted observed values 
are mostly conversions of the derived column densities with re¬ 
spect to Hj or CO column densities using a constant temperature, 
which introduces additional uncertainty. We discuss this aspect 
in detail in Sect. 14.21 


Garrod et al. (|2008|l suggested that cosmic-ray induced pho¬ 


todissociation of molecules (e.g. CH 3 OH) produce heavier radi¬ 
cals that effectively diffuse over grains at relatively higher tem¬ 
perature (>20 K) and produce other COMs such as C 2 H 5 OH. 
We found that a higher value of the cosmic-ray ionization rate of 
1.37 X 10 supports the formation of more COMs over grains 
and eventually also in the gas phase. At the initial stages of 
the evolution, the cr-16 models show nearly similar grain sur¬ 


face abundances for CjHgOH, CH 3 OCH 3 , 


and HCOOCH 3 over 


larger radii; but for the cr-17 models, these abundances not only 
differ by orders of magnitude but also vary with radii. The rela¬ 
tive abundance of these COMs with respect to CH 3 OH also dif¬ 
fer significantly between the cr-17 and cr-16 models (see Figs.|^ 
|A.lf|A.4| . During the warm-up phase, relative gas-phase abun¬ 
dances of COMs varies noticeably depending on the physical 
conditions. Moreover, even with a constant warm-up time scale, 
a slower rate of temperature increment provides favorable plat¬ 
forms to form more COMs (Tables 4] and |5]); a similar behav¬ 
ior was also found by [Garrod & Herbst ( [20()6[ ). In both models, 
the grain surface abundance of C 2 H 5 OH are higher than those 


of CH 3 OCH 3 


and HCOOCH 3 at the center, but the final gas 


phase abundance of CH 3 OCH 3 and HCOOCH 3 are higher than 
that of C 2 H 3 OH, suggesting that CH 3 OCH 3 and HCOOCH 3 also 
have effective gas-phase formation routes in addition to the grain 
surface formation routes. The grain surface abundance of these 
COMs reach their highest values at the outer regions where the 
highest temperature reaches only ~50 K. Interestingly, a similar 
pattern for grain surface abundances of HCOOCH 3 < C 2 H 5 OH 
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Fig. 4: Spatio-temporal abundance variation of selected COMs on the grain (dot-dashed line) and in gas phases (continuous line) 
are shown for the hot core model crl6-r3-\2. Black dashed and continuous lines show the radial variation of the density and the 
temperature respectively. Time-stamps (in year) for different temporal snapshots are annotated in the upper right corner of the 
respective plots. Note the alternative labeling of the density and temperature axes on the right side of the respective plots. 


Table 3: Abundance of COMs in Hot molecular cores: observed and simulated values 


Molecules 

IRAS 

16293-2422' 

Orion Hot 
Core^ 

NGC 

63341'' 

029.96^ 

G327.3-0.6'’ 

Garrod et al. 
(2008) 

Garrod 

(2013) 

CH 3 OH 

3.0(-7) 

2 . 2 (- 6 ) 

4.7(-6) 

6 . 6 (- 8 ) 

5.3(-8) 

3.2(-6) 

l.l(-5) 

C2HgOH 

- 

- 

4.7(-8) 

1 . 0 (- 8 ) 

2.6(-9) 

1^5(-9) 

5.9(-8) 

CH 3 OCH 3 

2.4(-7) 

6 . 8 (- 8 ) 

1 . 0 (- 6 ) 

3.3(-8) 

3.4(-7) 

3.0(-9) 

4.8(-8) 

HCOOCH 3 

2.0(-7) 

- 

2.4(-7) 

1.3(-8) 

5.0(-8) 

7.8(-ll) 

9.2(-8) 


References. (l) |Cazaux et al.|f2003| l; (2) |Crockett et al.|P014| l; (3) |Zernickel et al.]j2012^ ; (4) |Beuther et al.| j2009| l; (5) |Gibb et al.| ( |^00a| l;. 
Table 4: Simulated values of COMs at the center of various hot core models with the cosmic ray ionization rate of 1.37 x 10“^® 


Molecules 

r3-/2 

r3-/4 

i3-n 

r3.5-;2 

r3.5-;4 

r3.5-/7 

r4-/2 

r4-/4 

r4-/7 

r4.5-/2 

r4.5-/4 

r4.5-n 

CH 3 OH 

4.5(-7) 

1.3(-6) 

2 . 0 (- 6 ) 

1 . 0 (- 6 ) 

1.7(-6) 

1.7(-6) 

1.3(-6) 

9.8(-7) 

T9(-6) 

6.1(-7) 

1.3(-6) 

2.4(-6) 

C 2 H 5 OH 

4.3(-10) 

2.3(-9) 

9.5(-9) 

9.0(-10) 

3.8(-9) 

7.9(-9) 

Tl(-9) 

2.3(-9) 

8.3(-9) 

5.8(-10) 

3.0(-9) 

1 . 0 (- 8 ) 

CH 3 OCH 3 

6.6(-9) 

1.7(-8) 

3.5(-8) 

1.3(-8) 

2 . 0 (- 8 ) 

3.1(-8) 

1 . 6 (- 8 ) 

1.4(-8) 

3.2(-8) 

8.9(-9) 

1.7(-8) 

3.4(-8) 

HCOOCH 3 

3.4(-8) 

l.l(-7) 

1.4(-7) 

4.6(-8) 

l.l(-7) 

1.4(-7) 

4.7(-8) 

1.1(-7) 

T4(-7) 

3.7(-8) 

l.l(-7) 

1.4(-7) 
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Table 5; Simulated values of COMs at the center of various hot core models with the cosmic ray ionization rate of 1.37 x 10 


Molecules 

r3-/2 

r3-/4 

r3-/7 

r3.5-/2 

r3.5-/4 

r3.5-/7 

r4-Z2 

r4-/4 

r4-/7 

r4.5-/2 

r4.5-/4 

r4.5-/7 

CH 3 OH 

2 . 1 (- 6 ) 

2 . 2 (- 6 ) 

2.3(-6) 

2 . 2 (- 6 ) 

2 . 2 (- 6 ) 

2.3(-6) 

2 . 1 (- 6 ) 

2 . 2 (- 6 ) 

2.7(-6) 

2 . 1 (- 6 ) 

2 . 6 (- 6 ) 

2.3(-6) 

C,H 50 H 

9.3(-ll) 

1 . 2 (- 10 ) 

1.5(-10) 

9.3(-ll) 

1 . 2 (- 10 ) 

1.5(-10) 

9.2(-ll) 

1 . 2 (- 10 ) 

1.5(-10) 

9.4(-ll) 

1 . 2 (- 10 ) 

1.5(-10) 

CH 30 CH 3 

6.8(-9) 

7.0(-9) 

8.1(-9) 

6.4(-9) 

7.1(-9) 

8.1(-9) 

6.9(-9) 

7.0(-9) 

6.3(-9) 

6.9(-9) 

4.7(-9) 

8.1(-9) 

hcooch, 

3.4(-8) 

3.3(-8) 

3.7(-8) 

3.2(-8) 

3.3(-8) 

3.7(-8) 

3.5(-8) 

3.3(-8) 

2 . 6 (- 8 ) 

3.5(-8) 

1 . 8 (- 8 ) 

3.7(-8) 


< CH 3 OCH 3 < CH 3 OH has been observed at the outer regions 
of both the cr-16 and cr-17 models. During the warm-up phase, 
the gas phase abundance of CH 3 OCH 3 initially dominates other 
COMs; the final highest gas-phase CH 3 OCH 3 abundance ex¬ 
tends more than that of CH 3 OH in the cr-16 models, irrespective 
of density distribution and temperature evolution. In the cr-17 
models the spatial distribution of CH 3 OCH 3 shows two differ¬ 
ent peaks; the extent of the first peak is always narrower than 
the CH 3 OH distribution, but the extension and magnitude of the 
second peak depends on the density distribution (see Figs. |A.3| 
and |A.4| l. 


3.4. Spatio-temporal variation of formation routes of COMs 

As st ated earlier, we used the OSU gas-grain chemical net¬ 
work (|Garrod et al.||2008|l and the associated chemical path¬ 
ways are descnbed extensively in the corresponding paper 
(see also |Herbst & van Dishoeck| |2009| for a general sum¬ 
mary). For completeness, we brieily summarize the impor¬ 
tant grain surface reactions and outline the variation of forma¬ 
tion routes of COMs depending on the physical conditions. At 
the initial stages, the ice mantles primarily consist of s-H,0, 
S-CH 4 , S-H 2 CO, S-NH 3 , or S-CH 3 OH etc. S-H 2 O, S-NH 3 
and S-CH 4 are mai nly formed by successive hydrogenation of 
O, N, and C atoms (|Tielens & Hagen||1982|l and s-HtCO and 
s-CH„OH by successive hydrogenation ot CO (jCharnley &| 
Rodgers ||2009||. Cosmic-ray induced photodissociation of these 
molecules, specifically S-H 2 CO and S-CH 3 OH, produce other 
important radicals (chemical species that actively react with 
other species) such as s-OH, s-HCO, S-CH 3 , S-CH 3 O. During 
the warm-up phase (<100 K), these radicals effectively diffuse 
over grains and produce highly abundant oxygen-bearing COMs, 
for example, S-C 2 H 5 OH, S-HCOOCH 3 , and S-CH 3 OCH 3 ; the 
following reactions act as the main formation routes of these 
COMs depending on the physical conditions: 


s-H + S-CH3O — 1 S-CH3OH (1) 

s-H + S-HCOOCH3 —» S-CH3OH + s-HCO (2) 

s-OH -H S-CH3 —» S-CH3OH (3) 

s-HjCO -H s-CHjOH —» S-CH3OH -1^ s-HCO (4) 

S-CH2OH s-CO —» S-CH2OHCO 

S-CH2OHCO -H s-H —» S-CH2OHCHO ( 5 ) 

S-CH3 S-CH2OHCHO —» S-C2H5OH + s-HCO 

S-CH3 -I- S-CH2OH —» s-CjHjOH (6) 

S-CH3 -I- S-CH3O —» S-CH3OCH3 (7) 

S-CH3O -I- s-CO —1 S-CH3OCO 
s-H -H S-CH3OCO —1 S-HCOOCH3 

S-CH3O -I- S-H2CO —1 S-HCOOCH3 -I- S-H ( 9 ) 

S-CH3O -I- s-HCO —» S-HCOOCH3 (10) 


The mobility of the radicals strongly depend on the temper¬ 
ature, which shows a significant spatio-temporal variation and 
thus also affect the formation pathways of COMs. To facilitate 
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the discussion, we divided the hot core models into three dif¬ 
ferent subregions, central, middle, and outer regions. The typi¬ 
cal temperatures of these zones at the end of the evolution are 
>100 K, 60-100 K, and <60 K. In Table we summarize the 
important reactions at different regions that mainly produce the 
COMs described above. 

At the initial stages of evolution, that is at constant tempera¬ 
ture (15 K), Reaction[T]acts as the primary pathway of S-CH 3 OH 
formation for all the models. For the cr-16 models, up to 30 K, 
Reaction dominates and beyond that until the onset of evapo¬ 
ration (i.e., 100 K), Reaction [3] em erges as the main channel of 
S-CH 3 OH formation. Reaction!^ also occasionally becomes the 
main contributor of S-CH 3 OH formation in the middle region. 
However, in the outer region. Reaction [T] dominates up to 30 K, 
and beyond that. Reaction produces most of the S-CH 3 OH. 
These trends do not vary significantly with different luminosity 
evolution models. However, for the cr-17 models, some differ¬ 
ences are observed: although Reaction[T]initially dominates. Re¬ 
action [^begins to contribute soon and remains the main contrib¬ 
utory channel for temperatures up to 60-70 K. At higher temper¬ 
atures Reaction[^acts as the main formation route for S-CH 3 OH 
until it finally evaporates in the gas phase. But in the outer region 
Reaction [^begins to strongly dominate at relatively higher tem¬ 
peratures (> 15^0 K) also depending on the density. However, 
for temperatures around 60 K and above, the trend is similar to 
the cr-16 models that is Reaction [3 acts as the main formation 
pathway. 

Formation of other COMs, for instance, S-C 2 HJOH, 
S-CH 3 OCH 3 , and S-HCOOCH 3 , also strongly depend on the 
temperature since the diffusion of relatively heavy radicals such 
as S-CH 3 , S-CH 3 O, and S-CH 2 OH is involved. For the cr-16 
models the formation of S-C 2 H 3 OH is dominated by Reaction]^ 
around 20 K, and above that temperature until the desorption to 
the gas phase. Reaction [^becomes the main contributing chan¬ 
nel. For the cr-17 models. Reaction [^dominates until the mid¬ 
dle region up to 60 K and above that. Reaction acts as the 
main formation route; but in the outer region both of these reac¬ 
tions contribute to the S-C 2 H 3 OH formation. The formation of 
S-CH 3 OCH 3 is always dominated by Reaction ^ The formation 
of S-HCOOCH 3 in the cr-16 models is controlled by Reaction]^ 
and for temperatures lower than 60 K, and above that Reac- 
tionlUacts as the main formation route. For the cr-17 models the 
behavior is similar, except that Reaction does not contribute 
significantly to the formation of S-HCOOCH 3 . A more detailed 
discussion of the formation and destruction pathways of various 
molecular species including the COMs was reported in |Garrod| 
[etar] ( |2008l l. 

4. Analysis 

4.1. Qualitative comparison of simulated and observed 
spectra 

The most noticeable spectral features of high-mass star forming 
regions are the increase of the total number of spectral lines and 
their intensities with evolutionary stages (|Beuther et al.|[2009j 
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Table 6 : Spatial variation of formation routes of selected COMs 


Molecules 

Central region 

Radius:< 2000 AU 

Middle region 

2000^500 AU 

Outer region 
> 4500 AU 

S-CH 3 OH 

Reac. 1,2 &3 

Reac. 1,2, 3 &4 

Reac. 1 & 3 

S-C 2 H 5 OH 

Reac. 5 and 6 (>30 K) 

Reac. 5 and 6 (>25 K) 

Reac. 5 and 6 (>20 K) 

S-CH 3 OCH 3 

Reac. 7 

Reac. 7 

Reac. 7 

S-HCOOCH 3 

Reac. 8 , 10 & 9 (>60 K) 

Reac. 8 , 10 & 9 (>60 K) 

Reac. 8 & 10 



Frequency (GHz) 


Fig. 5: Simulated spectra of CH 3 OH for the hot core model crl6- 
r3-l2 at different evolutionary timescales. The adopted distance 
to the source and the telescope beam size is 2 kpc and 2 . 2 ". 



Fig. 6 : Simulated spectra of CH 3 OH for the hot core mod¬ 
els crl6-r3-\2, crl6-r3.5-\2, crl6-r4-\2, and crl6-r4.5-\2 at 
1.0 X 10^ year. The adopted distance to the source and the tele¬ 
scope beam size is 2 kpc and 2 . 2 ". 


Ceccarelli et al. 2010[ Gerner et al. 2014| l. Spectral variations 
are also observed in different sources with similar evolution¬ 
ary stages. Therefore, a characterization of the spectral evolu¬ 
tion based on the inline physio-chemical models can be used as 
a unique tool to constrain the evolutionary sequence of the ini¬ 
tial stages of high-mass star formation. Our simplistic model is 
unsuitable for any quantitative comparisons with the observed 
spectra, but it is expected that our model spectra at least repro¬ 
duce the general trends of spectral variation provided that our 
input hot core models resemble the conditions of the observa¬ 
tions. Comparing the simulated and observed spectra is thus also 
useful to determine the figure of merit of the input models. Mo¬ 
tivated by this aim, we simulated the spectra for various hot core 
models at different evolutionary stages at various wavebands and 
present the simulated spectra in Figs.|5]^ 

Simulated spectra of CH 3 OH for the hot core model crl6-r3- 
12 at different evolutionary stages are shown in Fig. the spec¬ 
tral range was chosen because it is often used in SMA/ALMA 
observations. In the literature the term line forest is often used 
when numerous complex molecules emit within a given fre¬ 
quency range and thereby create a forest of lines that perhaps 
hide emission from weaker molecules. In this paper we refer 
to the densely spaced ro-vibrational spectral lines of the COMs 
over a few GHz bandwidth as line forest. These lines arise from 
various transitions spanning over a range of energy states that 
trace different layers within the hot core. The spectra shown 
in Fig. cover two different line forests of CH 3 OH with rela¬ 
tively l^er (-60-300 K at -338 GHz) to higher (-350-650 K 
at -337 GHz) excitation energies. It can be easily seen from 
Fig. 1 ^ that the number of spectral lines and their respective in¬ 
tensities vary with evolutionary timescales. Moreover, it is also 
evident from the spectra that temporal evolution of spectral fea¬ 
tures are not simply a scaled-up version of the previous time- 


steps but originated from the combined variation of spatial distri¬ 
bution of temperature and abundances with evolutionary stages 


(Fig. C.li. Simulated spectra of CH 3 OH around 337 GHz for 
various hot core models with different density distributions are 
shown in Fig.|^ To avoid the self absorption features in the simu¬ 
lated spectra the highest value of the CH 3 OH relative abundance 
was restricted to 1 x 10 “^ (see Fig{^ compared with the val¬ 
ues shown in FigJ^ This discrepancy also indicates that CH 3 OH 
is overproduced in our models compared with the other COMs. 
This trend is also seen in the recent chemical models for hot 
cores ( |Garrod||2013a[ ). One possibility to explain this might be 
that there are some missing reactions in the chemical network 
that are useful to form other complex molecules using CH 3 OH. 
Since CH 3 OH is one of the parent molecules of other COMs, 
it would be interesting to compare the predicted relative abun¬ 
dance of other COMs with respect to CH 3 OH with the obser¬ 
vations to improve our understanding of the chemical evolution 
in hot cores ( [Neill et al.|2014| l. However, this discrepancy also 
emphasizes the importance of simulated spectra over the chem¬ 
ical abundances that are generally quoted by the astrochemical 
models to provide a better constraint on physio-chemical model 
of hot cores. 
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Fig. 7; Simulated spectra of selected COMs for the hot core model crl6-r3-l7 at 9 x 10“* year (left) and 1 x 10^ year (right) at 
different wavebands. The adopted distance to the source and the telescope beam size is 2 kpc and 2.2". 
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Fig. 8: Simulated spectra of selected COMs for the hot core model crl6-r3-l7 at 9 x 10"^ year (left) and 1 x 10^ year (right) at 
Herschel-HIFI Band 1. The adopted distance to the source and the telescope beam size is 2 kpc and ~40". 


Article number, page 10 of 


23 











































































































Article number, page 11 of 23 



-5.0 -2.5 0.0 2.5 5.0 


Fig. 9; Integrated line intensity maps of the CH 3 OH (v,= 0) [7(0,7)-6(0,6)] transition at 338.1245 GHz for the hot core model crl6-r3-\2 at 7.0 x 10"^, 8.0 x 10“*, 9.0 x 10^ and 
1.0 X 10^ year (left to right). The adopted distance to the source and the telescope beam size is 2 kpc and 2.2". 



Fig. 10: Integrated line intensity maps of the CH 3 OH (v,= 1) [7(4,4)-6(4,3)] and [7(4,3)-6(4,2)] transitions at 337.68559 GHz for the hot core models crl6-r3-\2, crl6-r3.5-\2, 
crl6-r4-\2 and crl6-r4.5-\2 (left to right) at 1.0 x 10^ year. The adopted distance to the source and the telescope beam size is 2 kpc and 2.2". 


mdhury et. al. 
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Fig. 11: top: Radial profiles of the CH 3 OH (v,= 0) [7(0,7)- 
6(0,6)] integrated intensity maps at 338.1245 GHz for the hot 
core model crl6-r3-\2 at various timescales. Bottom: Radial pro¬ 
file of the CH 3 OH (v,= 1) [7(4,4)-6(4,3)] and [7(4,3)-6(4,2)] in¬ 
tegrated intensity maps at 337.68559 GHz for various hot core 
models. The adopted distance to the source and the telescope 
beam size is 2 kpc and 2 . 2 ". 



arcsec 


Fig. 12: top: CH 3 OH 241.8877 GHz (E„= 72.5 K), and bottom: 
CH 3 OCH 3 241.9465 GHz (E„= 81.1 K) integrated line inten¬ 
sity maps of the hot core model crl6-r3-l7 at 1 x 10^ year. The 
adopted distance to the source and the telescope beam size are 
2 kpc and 0.3". 


Simulated spectra of different COMs such as CH 3 OH, 
CjHjOH, CH 3 OCH 3 and HCOOCH 3 of the hot core model 
crl6-r3-l7 for different frequency bands, evolutionary stages 
and angular resolution are shown in Figs. These spectra 

also follow the same trend: the specific line intensities become 
stronger and the spectra become more densely populated at later 
evolutionary stages. The frequency range s hown in Fig.[ 8 |is also 
covered by the Herschel-HIFI line surveys ( [Zernickel et al.|20T^ 
[Crockett et al.|2014| l. A qualitative comparison of the simulated 
and observed spectra reveals that the relative line intensities of 
CH 3 OCH 3 and HCOOCH 3 with respect to CH 3 OH are compa¬ 
rable with the observations (Fig. EIJ- 

Integrated line intensity maps of the CH 3 OH transition at 
338.1245 GHz at different timescales are shown in Fig. [9] The 
corresponding radial profiles for the integrated intensity maps 
are shown in Fig. (top panel). In both figures, the radial pro¬ 
files of molecular emission expand with time. This is due to the 
radial expansion of the gas phase abundance of CH 3 OH with 
evolutionary time resulting from the increase in the stellar lumi¬ 
nosity (Fig.[^a)) and the successively higher mass that is heated 
to higher temperatures. However, as our models use a static den¬ 
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sity distribution and do not yet include dynamics, the spatial ex¬ 
tension of temperature, abundance, emission maps, etc. should 
be considered as approximate values. In Fig.[T^ integrated line 
intensity maps of CH 3 OH transition at 337.68559 GHz for vari¬ 
ous hot core models are shown. The corresponding radial profiles 
are shown in Fig. (bottom panel). The variation in the emis¬ 
sion maps can be attributed to the different density distributions 
of these models since the other input parameters are the same for 
all the models. These maps demonstrate the importance of the 
underlying physical stracture on chemical evolution and spectra 
of hot molecular cores. It also vouches for the potential of this 
study to achieve a better understanding of the hot cores and thus 
the high-mass star formation scenario. 

One of the important parameters of astrochemical models is 
the desorption energy of the respective molecule^ it also con¬ 
trols the gas-phase abundance of COMs. Recent observational 
studies also started to compare the spatial distribution of COMs 


® A list of binding energies for various molecules can be found 
at http://www.astro.comell.edu/~rgarrod/wp-content/uploads/2013/03/ 
GH06_binding.txt 
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Fig. 13: Comparison of the input abundance distribution of the hot core model crl6-r3-l7 and the single-point physical parameters 
obtained by the myXCLASS analysis. 
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using high resolution emission line maps ( |Oberg et al. 2013| l. 
The relative gas phase abundance of CH 3 OCH 3 is more extended 
than that of CH 3 OH in the hot core models, therefore, it will be 
interesting to verify whether this trend is also seen in observa¬ 
tions. We show the integrated intensity maps of the CH 3 OH: 
241.8877 GHz (E„= 72.5 K) and CH 3 OCH 3 : 241.9465 GHz 
(E„= 81.1 K) transitions for the model hot core crl6-r3-l7 at 
1 X 10^ year in Fig. 12 This hgure demonstrates that the bond 
strength of the molecule to the grain can influence the size 
of the emissive region (i.e., a lower desorption energy places 
the molecule in the gas phase over a larger surface area). We 
can compare the observed emission maps from a range of par¬ 
ent molecular species, each with different desorption energies, 
but sampling transitions with comparable excitation conditions. 
Then these emission prohles can be used to constrain the rela¬ 
tive desorption energies of the associated molecules. Moreover, 
with grounding to laboratory work for a particular molecule, we 
might be able to place this on an absolute scale. 

It has been found that the abundances of various COMs 


Table 7: Physical parameters of the hot core model crl6-r3-l7: 
myXCLASS analysis 


Time Molecule 


myXCLASS analysis 

(year) 

Frequency 

(GHz) 

Source Tempera- Relative 

size (") -ture(K) abundance 


CH3OH 

248-251 

2.6 

139 

1.2 X 10-“ 


337-339 

1.86 

137 

1.2 X 10-“ 

C,H,OH 


1.42 

139 

2.2 X 10-* 

CH3OCH, 

332-336 

1.75 

126 

6.3 X 10-* 

HCOOCH3 


1.16 

194 

1.4 X 10-* 

CH3OH 

248-251 

>2.2 

144 

1.8 X 10-“ 


337-339 

>2.2 

141 

1.9 X 10-“ 

C,H,OH 


>2.2 

133 

4.5 X 10-* 

CH3OCH3 

332-336 

1.85 

166 

8.5 X 10-* 

HCbOCH, 


1.75 

189 

7.5 X 10-* 


dances using H 2 or CO column densities at constant temperature. 
We used the total hydrogen density to estimate the relative abun¬ 
dance of the respective molecules (see Tables and [^. There 


in crl7 models are typically lower than that of crl6 models, 
and these values are unsuitable to reproduce the observed line 
strength of different COMs. Since our current approach is based 
on static density distribution, it is not possible to distinguish the 
effect of dynamics vs. cosmic-ray ionization rates on the ini¬ 
tial ice formation and hence on the hnal gas phase abundance 
of COMs. More detailed studies are needed to constrain the 
cosmic-ray ionization rate at the sites of high-mass star forma¬ 
tion. 


4.2. myXCLASS fit to the simulated spectra 


Many of the line-survey studies used XCLASS (SchiUce et al. 


is another notable difference between the XCLASS analysis and 
this work: the hot core models predict the temperature and abun¬ 
dance prohles for hot cores, whereas XCLASS and other meth¬ 
ods can only derive beam-averaged values. Therefore, it would 
be interesting to compare these single-point physical parameters 
with the spatial distribution predicted by our hot core models. 
We carried out an analysis to re-estimate the physical parame¬ 
ters from the simulated spectra using the myXCLASS interfac^ 
(Moller et al. 2014, in preparation) for CASAp^ the interface 
also includes the model optimizer package MAGIX (Modeling 
and Analysis Generic Interface for external numerical codes) 

(Moller etal pon I 

^Fl., 


We htted the simulated spectra of CH 3 OH, 


.. 2 H 5 UH, CH 3 UCH 3 , and HCOOCH 3 convolved with a tele- 


IW 

20011 IZernickel et al.||2012l |Crockett et al.||2014l |Neill 

® https://www.astro.uni-koeln.de/projects/schilke/ 

et al. 

2014|l or similar programs to estimate the physical param- 

myXCLASSInterface 


eters from the observed spectra and derived the relative abun- http://casa.nrao.edu/ 
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scope beam size of 2 . 2 "at two different evolutionary stages over 
selected frequency ranges. The fitted spectra are overplotted on 
the simulated spectra in Fig. |E.1| We list the physical param¬ 
eters obtained from the myXCLASS analysis in Table. |7j the 
relative abundances were derived using the hydrogen column 
density, 1.05 x 10^"^ cm“^, as estimated for the hot core models 
hmc-r3. We also added white noise to our simulated spectra and 
htted the noisy spectra with myXCLASS. We do not hnd any 
significant difference in the htted parameters, indicating that the 
myXCLASS htting values are robust against noise. We also plot 
the relative abundances as a function of temperature for differ¬ 
ent molecules at different evolutionary timescales and overplot 
the single-point physical parameters estimated by myXCLASS 
analysis in Lig. 13 The myXCLASS estimates are compara¬ 


ble with the physical parameters at the jump radius of these 
molecules (where the gaseous abundances sharply decrease) (see 
also Lig. |B.l| i; these regions have the largest area-hlling fac¬ 
tor in the beam, therefore, the results are expected. However, 
myXCLASS underestimates the temperature of the hot cores. 
Our analysis indicates that the typical temperatures of the hot 
cores are a few times higher than the canonical value of 100 K, 
which is also consistent with recent observations ( |Sridharan et al.| 
|2014| l. These results also suggest that a comparison of the simu¬ 
lated spectra of the 3D hot core models with the high resolution 
observations (e.g., ALMA) will reveal a more realistic physio- 
chemical structure of the hot cores than the simple htting meth¬ 
ods. 


5. Summary and discussion 

Molecular line observations are important probes to reveal the 
chemical reservoir of hot cores ( [Schilke et al.|[T9971 | 2001 | l. 
These observations are also useful to explore the physical condi- 
tions and chemical evolution of high-mass star forming regions 


( Beuther et al.||2069 Zernickel et al.|[2012 Oberg et'Sl]|2013 i. 
Recent studies also used single dish observations of a large num¬ 
ber of sources at various evolutionary stages to investigate the 
chemical evolution of high-mass star forming regions ( |Hoq et ah] 
2013[ Gerner et al.|2014 1. Selected studies also explored the ra¬ 
dial distrihution of density and temperamre along with chemi¬ 
cal evolutionary models to understand the underlying physical 
structure, abundance evolution , evolutionary stage s, etc. (Doty 
et al.||20(H Oberg et al.||2013 Gemer et al.||2014 i. In general, 
these studies analyzed the observations by assuming a density 
distribution guided by empirical studies and an estimate of the 
thermal structure that would result from the density profile and 
a source with a particular luminosity. Doty et aT] (|2006)l deter¬ 


mined the temperature evolution including the luminosity evolu¬ 
tion for high-mass protostars, but did not include any grain sur¬ 
face diffusion reactions, which are cracial for COM formation 
on grain surface. 

In this work we have used different temporal evolution pro- 
hles of the luminosity for a forming high-mass star. The den¬ 
sity distribution of the hot core models were setup using recent 
dust-continuum observations assuming a gas to dust mass ra¬ 
tio of 100. Within this framework, the thermal evolution was 
determined with the radiative transfer code RADMC3D based 
on the adopted protostellar luminosity evolution. This then re¬ 
sulted in a time-variable thermal prohle that has a profound ef¬ 
fect on the chemistry through sublimation and grain surface dif¬ 
fusion, subsequent gas phase processing, etc. We used the chem¬ 
ical evolutionary code Saptarsy along with a gas-grain chemi¬ 
cal network to explore the spatio-temporal evolution of molecu¬ 
lar abundances in hot cores. We then explored the evolution of 


abundance distribution of selected oxygen-bearing COMs such 
as CH 3 OH, CjHjOH, HCOOCH 3 , and CH 3 OCH 3 in detail un¬ 
der varying physical conditions. We showed that the spatial dis¬ 
tribution of gas phase abundances of the COMs expand with the 
increasing temperature of hot cores. These prohles provide a de¬ 
tailed variation in abundances against the physical structure of 
the hot cores. The average behavior of these prohles is similar to 
the so-called jump models (two different abundances are used for 
the upper and lower temperature regime relative to a characteris¬ 
tic temperature, typically 100 K) for species that primarily form 
on the grain surface and are transformed to the gas phase mainly 
via thermal desorption. In the evolutionary model the jump ra¬ 
dius or the evaporation font (where the temperature is equal to 
the characteristic temperature of the jump abundance) also ex¬ 
pands with time. We also showed the variation of the evapora¬ 
tion front associated with different temperatures for the hot core 
model hmc-r3-l7 (see Fig. |B.l| i, which will be useful to set the 
jump abundances for any particular evolutionary stages. More¬ 
over, the predicted abundance profiles can be used as templates 
to model the observed spectra of hot cores, but one should also 
keep in mind that these prohles vary signihcantly depending on 
the adopted physical structure and luminosity evolution. 

Emission lines of COMs (often referred to as weeds) are 
one of the salient observational features of hot cores. These 
molecules are highly abundant (> 10 “^ with respect to hydrogen) 
and have rich ro-vibrational emission spectra. A signihcant im¬ 
provement over previous studies was achieved by simulating the 
spectra for hot core models that can be directly compared with 
observations. The simulated spectra for different hot core models 
were generated by using RADMC-3D at different spectral wave¬ 
bands to explore the spectral evolution in hot cores. The qualita¬ 
tive comparison of the observed and simulated spectra showed 
that the hot core models successfully reproduce the observed 
trends, that is the increase of the total number of emission lines 
and their associated intensities with evolutionary timescales. The 
adopted timescales for the chemical simulation of the cold core 
phase (~ 10 ^ year) or the hot core lifetime (~ 10 ^ year) are con¬ 
sistent with recent independent studies (Tackenberg et al.|2012 


[Gerner et al.|2014| l. These results indicate that the self-consistent 
models present a reasonable evolutionary scenario for hot cores. 
Based on these evolutionary models, we discuss some additional 
features that can be compared with observations. 

The hot core models predict that the emission line maps of 
various transitions of different COMs also expand with time. 
With increasing temperature, the material around the central pro¬ 
tostars heats up, and consequently, the gas phase abundance of 
various COMs also spreads out, depending on the binding en¬ 
ergies of the respective molecules. These processes collectively 
contribute to the expansion of the emitting regions of various 
transitions. We argue that for a given molecule the line forest 
(see Sec. |4. l| l provides detailed information about the excitation 
conditions, which can be used to extract key information regard¬ 
ing the binding energies of the molecules to the grain surface, 
abundance prohles, etc. This will provide further constraint to 
the chemical models that sometimes suffer from lack of labora¬ 
tory measurements of the binding energies. A satisfactory match 
between the observed and simulated spectra along with the emis¬ 
sion line maps thus provides a well-constrained physio-chemical 
structure of the hot cores, which is a signihcant improvement 
over existing htting methods that estimate the beam-averaged 
single-point physical parameters. High resolution observations 
(such as from ALMA) for a statistically signihcant samples of 
high-mass prestellar cores will be available soon. The tempera¬ 
ture structure of these cores can be re-constructed by using ra- 
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diative transfer modeling and can be further constrained by com¬ 
paring the simulated emission line maps with the observations. 
The fitting process also place some constrain on the luminosity 
of the central stars, gas phase abundance profile and the desorp¬ 
tion energy of the associated molecules. These models will be 
particularly useful for the molecules that primarily form on the 
grain surface and are transformed to gas phase mainly by ther¬ 
mal desorption. Moreover, modeling the high resolution obser¬ 
vations of the high-mass protostellar cores with similar masses 
but at different evolutionary stages will be useful to estimate the 
protostellar luminosity function for high-mass protostars, which 
is an input parameter of these models. In summary, the full po¬ 
tential of these models can be explored by using the observations 
that can probe scales of few thousands AU. Currently, it is also 
possible to estimate the evolutionary stages using these models, 
but the estimated age should be considered with caution because 
a static density distribution is used. With a temporal variation of 
density distribution, it may be possible to diagnose the spectral 
features or abundance ratios that can be used as a proxy for the 
evolutionary stages, similar to the so-called chemical clocks. 

Along with the variation in chemical abundances, the dif¬ 
ferences in formation and destruction pathways of various 
molecules as a function of physical conditions can be explored 
using our model. This is specifically helpful to verify the pro¬ 
posed chemical evolutionary scenarios. CH^OH acts as the par¬ 
ent molecule for the other COMs that are considered in this 
work. The spatio-temporal variation of the abundance ratio of 
these COMs with CH 3 OH can be used to verify the reaction 
pathways discussed in Sect. 3.4 The abundance and spectra of 
intermediate species such as CH 3 O and CH 2 OH will also be use¬ 
ful to constrain the formation scenario of COMs. 


6. Conclusion 

We investigated the formation and evolution of COMs in HMCs 
by constructing 3D physio-chemical models guided by recent 
empirical and modeling studies of high-mass star formation. We 
combined radiative transfer calculation with chemical evolution 
to explore the spatio-temporal evolution of the COMs in hot 
cores and also to generate the synthetic spectra at selective wave¬ 
bands. The grain surface abundances of COMs are converted 
into gas phase, where the temperature is higher than the water 
ice evaporation temperature (100 K). The gas-phase abundance 
profile of various COMs will furthermore be useful for modeling 
the spectra of individual sources. We also find that jump abun¬ 
dance profiles are still a good approximation for the COMs that 
mainly form on the grain surface and return to the gas phase by 
thermal desorption (e.g., C 2 H 3 OH). We predicted the temporal 
variation of the jump radius that can be used as a template to set 
the jump abundances according to the evolutionary stages. We 
successfully reproduced the trends in spectral variation observed 
in hot cores within the typical lifetime of HMCs of 1 x 10^ year. 
These results indicate that the adopted physical parameters and 
the chemical model ( |Garrod et al.||2008l l are realistic estimates 
of the observed hot cores. We also discussed that a satisfactory 
match between the simulated and observed spectra obtained by 
iterating over the input parameter space will be useful to deci¬ 
pher the physio-chemical structure of the hot cores. This tech¬ 
nique when applied to model the observation of a handful of hot 
cores can also constrain the luminosity evolution of the central 
protostars and the binding energies of the molecules to the grain 
surface. Our approach is more effective in extracting the phys¬ 
ical parameters (such as temperature) using observations than 
the simple spectral fitting methods that estimate single-point val¬ 


ues for the whole source. The modeling approach discussed in 
this paper will be particularly helpful to reconstruct the physical 
structures of hot cores and to constrain the timescales and chem¬ 
ical evolution during the initial stages of high-mass star forma¬ 
tion. 
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Appendix A: Spatio-temporal variation of selected 
COMs in various hot core models 

Appendix B: Temporal variation of the radii with 
temperature above 100 K, 150 K, and 200 K 

Appendix C: Comparison of simulated spectra for 
different evolutionary timescales 

Appendix D: Comparison of simulated and 
Herschel-HIFI spectra 

Appendix E: myXCLASS fit to the simulated spectra 
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Fig. A.l: Same as Fig. but for hot core models crl6-r3-l7 (upper panel) and crl6-r3.5-l4 (lower panel). Note the alternative 
labeling of the density and temperature axes on the right side of the respective plots. 
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Fig. A.2; Same as Fig. but for hot core models crl6-r4-l2 (upper panel) and crl6-r4.5-l4 (lower panel). Note the alternative 
labeling of the density and temperature axes on the right side of the respective plots. 
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Fig. A.3: Same as Fig. but for hot core models crl7-r3-l2 (upper panel) and crl7-r3.5-l4 (lower panel). Note the alternative 
labeling of the density and temperature axes on the right side of the respective plots. 
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Fig. A.4: Same as Fig. but for hot core models crl7-r4-l7 (upper panel) and crl7-r4.5-l4 (lower panel). Note the alternative 
labeling of the density and temperature axes on the right side of the respective plots. 
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Fig. B.l; Temporal variation of the radii with temperature above 100 K, 150 K, and 200 K. Radii of the spatial distribution of various 
molecules as derived from the myXCLASS analysis (see Table|7]l are overplotted with different colors. The arrow indicates that the 
radii of CH 3 OH and CjHjOH at 1 x 10^ year are only lower limits. 



Fig. C.l; Simulated spectra for the hot core model crl6-r3-l7 at 9 x 10^ year, overplotted on the spectra at 1 x 10^ year. The 
intensity of various spectral lines at two different epochs cannot be matched using a scaling factor. Spectral changes originate from 
the variation of temperature and abundances at different timescales. 



Fig. D.l; Scaled-up version (multiplied by a factor of 2) of simulated spectra for the hot core model crl6-r3-l7 at 9 x 10^ year, 
overplotted on the Herschel-HIFl spectra of NGC63341. 
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Fig. E.l: Fitted spectra obtained from the myXCLASS analysis (solid black line), overplotted on the simulated spectra (colored 
dashed lines) of different COMs. 































